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Abstract. Using the decay of the out equilibrium spin-spin correlation function we 
compute the equilibrium Edward- Anderson order parameter in the three dimensional 
binary Ising spin glass in the spin glass phase. We have checked that the Edward- 
Anderson order parameter computed from out of equilibrium numerical simulations 
follows with good precision the critical law as determined in experiments and in 
numerical studies at equilibrium (which allow us to estimate the /3 critical exponent). 
Finally we present a large time study of the off-equilibrium fluctuation-dissipation 
relations and we find strong discrepancies (in the low temperature region) between the 
numerical data and the droplet theory predictions and agreement with the predictions 
of the replica symmetry breaking theory. 
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1. Introduction 

The characterization, using numerical simulations, of the phase transition in the three 
dimensional Ising spin glasses has been a challenging problem. Recently a clear picture 
of the phase transition and good estimates of the critical exponents have been obtained 
for both Gaussian and bimodal disorder by working at equilibrium [U El E! • 

However a characterization of the phase transition using out of equilibrium 
techniques is still lacking (see reference jl] for a detailed discussion). In the first part of 
this paper we will address this problem (simulating the bimodal disorder). In particular 
we will compute the order parameter using out of equilibrium techniques pi and we 
will characterize the transition using this observable. In addition we will confront our 
data with previous estimates of the critical point and critical exponents for this model 
(obtained from numerical simulations and from experiments). The behavior of this 
observable will permit us to discard (again) a Kosterlitz-Thouless like phase transition 
(as done in equilibrium jl] , that we will refer in the following as XF-like scenario) for the 
transition [l]. Moreover, we have studied the dependence of the order parameter with 
the size of the system. Hence, we will present on this paper the first direct numerical 
computation of the Edwards- Anderson order parameter in the three dimensional Ising 
spin glass (obtained out of equilibrium). 

This kind of study was performed in the past in four dimensions J2j (see also [3 IE!) 
but is still lacking in three dimensions (the interesting physical dimensions). 

The second part of the paper is devoted to the study of the fluctuation-dissipation 
theorem out of equilibrium. This kind of analysis have attracted a large amount of work 
(analytical, numerical and experimental) in the last years IIUI ITTl IT^ IT^ I14j . 

Using the results of reference jHj and assuming that the three dimensional Ising 
spin glass presents stochastic stability (until now it has not been rigorously proved 
but there are numerical evidences ^^I) one can relate the fluctuation-dissipation curves 
with equilibrium properties and so, compute or measure the equilibrium probability 
distribution of the overlap. This computation or measurement is very important since 
it should discern between the different theoretical approaches in competition, which try 
to describe the behavior of finite dimensional spin glasses (e.g. the Replica Symmetry 
Breaking (RSB) approach [TT)| IT7j or the droplet model[TH]). 

The goal of this (last) part of the paper is twofold. First, to check if the order 
parameter computed in the first part of this paper matches well in the fluctuation- 
dissipation (FD) curves. This is important since this value marks the point in which the 
FD curve departs from its pseudo-equilibrium regime, and the behavior of the curve from 
this departing point is a clear fingerprint whether or not the system behaves following 
the RSB theory or the droplet model. 

And the second goal is to study the finite time behavior (for really large times) 
of the curves in order to see how the asymptotic form of the FD curves is built up. 
This is important, since until now, the numerical simulations ^2j and experiments 
[T3] show up a behavior compatible with the Replica Symmetry Breaking description 
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and incompatible with droplet theory. One can argue that the curves reported in the 
literature ^3 ^] are not asymptotic and that the asymptotic curve is compatible with 
droplet theory and no compatible with RSB. 
Finally, we will report the conclusions. 



2. The model and Numerical simulations 



We have simulated a three dimensional system in a cubic lattice with helicoidal boundary 
conditions of size L and volume V = L^. The Hamiltonian is 

'H = - Jij(^i^j ' (1) 

<«j> 

where < z,j > denotes the sum over the first nearest neighbors, (jj = ihl are Ising 
variables and = ±1 are quenched random variables with a bimodal probability 
distribution with zero mean and unit variance. We have used the standard heat-bath 
algorithm (local dynamics) to simulate the three-dimensional lattice. 

We will introduce the observables measured in our work. Firstly, the order 
parameter (the Edwards Anderson one) is defined as: 



gEA=(fT.)2, (2) 



where, as usual, we use ((■ ■ ■)) and (■ ■ ■) to denote thermal and quenched disorder 
average respectively. 

In addition, the spin-spin correlation function has been computed using 

1 ^ 

C{t,t^) = -Y,a^{t)<Ji{t^) . (3) 
^ 1=1 

We can obtain formally the order parameter from this correlation as the double limit: 
gEA = lim lim C(t, t^,) . (4) 

t — >00 t-w — ^OO 

Notice that the order of the limit is crucial in obtaining the order parameter. We will 
use this equation to extract gsA from the out-of-equilibrium data. 

We will study in the last part of the paper the finite time behavior of the violation 
of the fluctuation-dissipation relation in the three dimensional spin glass. We will review 
shortly the main equation of the off-equilibrium fluctuation-dissipation equations (see 

for more details): 

R{t,M) = \x{C{t,,t,)f-^^^ , (5) 

where, ti > t2, -R(ti, ^2) is the response of the system to the magnetic field perturbation 
(i.e. the magnetic susceptibility of the system: R{ti,t2) = m(ti,t2)/h) and X{C) is 
the, in principle unknown, function which controls the violation of the fluctuation- 
dissipation theorem. Integrating this equation in ^2 and taking the perturbing field as 
h{t) = h6{t — tuj) we finally obtain (working in the linear-response region): 

m{t) ~ ph f du X{u) . (6) 
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In the regime ti ^ t2 ^ 1 we reach the equihbrium, and it is possible to show 
that C(ti,t2) — > q. In addition X{q) — > x{q) = J^^^.^dq'P{q'), where x{q) is the integral 
of the probability distribution of the overlap at equilibrium fH]- Hence, in this regime 

cni HI ca 1131111, 

m(t) c:^ 3h du x(u) . (7) 

Jc{tu) 

Furthermore, we can define 

S{C) ^ [' dq x{q) , (8) 



SiC{t,tJ). (9) 



so, 

m{t)T 
h 

Both, in droplet theory and RSB (see reference j2l], in particular its figure 10), S{C) 
is the straight line 1 — C for C G [^ea, !]• However, for C < ^ea the behavior is very 
different: in the droplet theory S{C) is constant in this region and in RSB S{C) is a 
growing function with curvature. We recall that knowing the initial point, S{C = 0), 
we can compute ^ea in the droplet theory as 

^droplet ^ ^ _ ^ _ (^0) 

This technique allows us to compute, taking the appropriate limit, the equilibrium 
function x{q). 

Finally, we report that all the numerical simulations have been obtained with the 
SUE machine This is a dedicated machine, designed for the simulation of the three 
dimensional Edwards- Anderson model with first neighbour couplings |16j. the system 
that is being studied in the present work. It consist of 12 identical boards. Each single 
board is able to simulate 8 different systems, updating all of them at each clock cycle. 
SUE reaches an update speed of 217 ps/spin with a clock frequency of 48 MHz. The 
on-board reprogrammability permits to change in an easy way the lattice size, or even 
the update algorithm or the Hamiltonian. The SUE machine is connected to a Host 
Computer running under Linux. SUE is in charge of the update of the configurations, 
and the host computer is in charge of measurements and analysis. The main electronic 
devices of each SUE board are the Altera family, that performs the update. Other 
devices store the spins and couplings variables. One of the Alteras is devoted to generate 
random numbers in a fast way (for more details, see Ref. ^5). Up our knowledge, 
SUE has been the fastest dedicated machine in the simulation of the three dimensional 
Edwards- Anderson model. 



3. Computation of the Edward- Anderson Order Parameter 

In order to compute the Edward- Anderson order parameter (^ea), we have carried out 
several runs for two lattice sizes and different temperatures: (3 = 1/T = 2.00, 1.67, 1.25, 
1.05, 1.00, 0.95 and 0.91 for L = 30; and p = 2.00, 1.67, 1.25 and 1.00 for L = 60 . For 
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Figure 1. Out of equilibrium spin-spin correlation function C{t,tw) computed for 
i = 60 and = 1.25. Top: C{t,t^) versus time, t. Bottom: C{t,t^) versus waiting 
time, i^, obtained by studying figure in top for several fixed times t in order to find 
the limit — > oo behavior of C{t, t^). The continuous lines in the plot arc the fits to 
equation Ullfl . Notice that for the curves with larger waiting time we have chosen to 
show not all the fits to Ullfl in order to present a clean figure (the quality of the fits is 
the same for all the waiting times). 

all of them we have averaged over 58 samples. In figure (Q) we report the curves C(t, t^) 
as a function of time t. 

We have checked that the behavior of C(t, t^) for ^ 1 follows with high precision 
the behavior (as in higher dimensions, see [HIIZ!; this is just an Ansatz): 

C(t,U = a(t) + 6(t)CW, (11) 

where a{t) is related with the value of gsA- In order to find it out we have first obtained, 
from figure ((T)) top, the curves C(t,tw) vs. t^ for several fixed values of t (typically, 
from 8192 to ~ 3.7 x 10^ Monte Carlo steps) (see figure dH) bottom). We have fitted 
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Figure 2. Function a{t) defined in equation 1)11(1 computed for L = 60 lattice size and 
/3 — 1.25 (top) and /3 = 1.67 (bottom). Notice that the last points of the curve can 
be fitted to a constant in the two plots. In addition we have drawn the function a{t) 
with only 16 samples in order to show the dependence of the extrapolated value (i.e. 
the plateau) on the number of samples. 



these curves to the functional form defined in (fTT|) obtaining in this way the behavior 
of a{t) as function of t (we show these fits in figure IQ) . From a{t) and for t ^ 1, we 
can obtain the value of qea (since asymptotically a{t) must became ^ea)- To achieve 
this aim, we have fitted the last points of a{t) versus t to a constant function (since a{t) 
shows a clear plateau, see Fig.(j21)). In this way, we have implemented the double limit 
in equation (0}. The results obtained from these fits are shown in Fig. ((21) • 

We have checked that for (3 > 1.00 the values for ^ea are the same for both L = 30 
and L = 60. In (3 = 1.00 the difference is about 1.5 standard deviations. In addition we 
have run a L = 20 lattice at /? = 0.91 and f3 = 1.00: these data show finite size effects 
as expected since they lie near the critical point (see figureij^I)). 
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4. Characterizing the Phase Transition 

As we mentioned before, we have checked that the ^ea, which we have computed out of 
equihbrium, follows with good precision the critical law of the order parameter 

qMf^) = A{(3 - f3,r\ (12) 

where we have denoted (3q the usual (3 exponent of the order parameter (in order to 
avoid confusion with the usual notation (3 = 1/T) 

By fitting only the points closer to the critical one (satisfying (3 < 1.25) we obtain 

pc = 0.866(2) = 0.52(9) , (13) 

with a x^/d.o.f = 1.13. This figures compare really well with the numerical values 
obtained at equilibrium P, namely: Pc = 0.88(1) and Pg = 0.71(5). In particular the 
difference between the two estimates of Pg is 0.19(11), less than two standard deviations. $ 

In addition, we can compare with experiments. In reference [SHI was found 
Pg = 0.54(10)11 which is in a very good agreement with our out equilibrium value. 




n I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I u 

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 

P 

Figure 3. q'^^ versus (3 for three lattices sizes L = 20,30 and 60. Tlie continuous 
line is the fit reported in the text. 

X Notice that in reference ^ corrections to scaling were taken into account. In our estimate there is 
no scaling corrections, hence our error are smaller than the error quoted in 1 : i.e. our error bars are 
underestimated. 

§ See also for a non Universality scenario: they reported /3c — 0.84(1). 
II Note that both results in ^ and |2S1 come from different methods. 
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We have also checked that ^ea follows with good precision the critical law 

{qEAm'/^^ = A{(3-(3,), (14) 

Again, we have only used in the fit the points with (3 < 1.25 (critical region). Moreover 
we can fix f3q to the experimental value, obtaining again a compatible value with the 
equilibrium one: Pc = 0.8603 (6) (236), where the first error is statistical and the second 
error comes from the error of the experimental jSq. In addition, by fixing jSq to the 
numerical simulations value we obtain jSc = 0.820(3)(13), less than three standard 
deviations from the numerical value. 

All figures reported in this analysis are compatible with latest estimates of the 
critical exponents. In reference [2] /5c = 0.893(3) and jSg = 0.723(25) were reported. 
In addition a diluted version of this model was studied in jH] and Pq = 0.723(50) was 
reported. 

Finally, we remark that our numerical results from both f3c and f3q must suffer from 
the systematic error coming from the dependence of qea with L near the critical point 
(as shown the L = 20 runs). At (3 = 1.00 we have three different values of the order 
parameter that fit to the law 

gEA(^) = gEA(oo) + , 

where b and c are constants. This is the finite volume correction equation which holds 
in the low temperature phase ^. We have obtained c = 3.54 and gEA(oo) = 0.49 (notice 
that we are fitting three points to a three parameter function) to be compared with 
gEA(^ = 60) = 0.485(6) and gEA(^ = 30) = 0.47(1). At (3 = 0.91 (the nearest value we 
have to the critical point) we have only two points, that anyhow, we can try to fit to 
equation (jH) fixing c = 3.54, obtaining gEA(oo) = 0.278 (no error bars can be reported 
since, again, the number of degrees of freedom in this fit is zero) to be compared with 
the value of our largest lattice Qea^L = 60) = 0.26(1), so this limited analysis suggests 
that the L = 30 lattice is asymptotic in its error bars in the region /3 > 0.91. Hence, we 
are confident that our final estimates of (3c and (3q should have small systematic error 
coming from finite size effects. 

We remark that testing the dependence of ^ea with the lattice size, for large lattices 
(e.g. L = 60) near the transition is not accessible even using the SUE machine. 

5. Finite Time Effects in the Fluctuation-Dissipation relations 

We have performed several runs again with SUE machine, in a lattice of size L = 60 for 
different temperatures: /? = 1.25, 1.10, 1.05, 1.00 and 0.95. We have used the following 

*li In reference was checked that in the three dimensional Gaussian Ising spin glass the position of 
the maximum of the equilibrium probability distribution of the overlap follows this law with c = 1-5(4) 
by fitting L < 16. Notice that in our case we are using 20 < L < 60 data and we simulate the ±J 
model and that the c exponent could depend on the temperature. Notice that usually in equilibrium 
small lattices develop larger order parameter, however, in our dynamical approach we have found the 
opposite behavior. 




Figure 4. Fluctuation-dissipation curve out of equilibrium for one of the lowest 
temperature simulated /3 — 1.25, L 60 and one waiting time for two different 
perturbing magnetic fields, h — 0.01 and 0.03, in order to check linear-response. 



standard procedure. We let the system evolve during a time t^, just after this time, 
a field h = 0.03 is plugged, seeing the response of the system and recording the 
magnetization and the correlation function. Then it is possible to extract the value 
of Qea, for the particular f3 being analyzing at that moment, from the point where 
the curve leaves the linear regime, that is, where mT/h does not follow the pseudo- 
equilibrium line (1 — C)/T. 

The choice of the field strength applied to the system has not been arbitrary. We 
need to stay in the linear-response region. We have checked this by simulating different 
magnetic fields: h = 0.01,0.03,0.05 and 0.10. Finally we have selected a safe value for 
h: h = 0.03, which is a compromise between large and small fields (notice that small 
fields induce strong noise in the measures). In figure (jH) we have shown the FDT curve 
for a waiting time and two perturbing magnetic fields {h = 0.01 and 0.003) in order to 
test that we are in the region in which linear-response holds. It is clear from this figure 
that the curve, inside the error bars, is independent of the perturbing magnetic field. 

In the droplet model, the curve X{C) departs horizontally from the straight line 
1 — C, the final value of the horizontal fine being mg^y^T/h (i.e. S{C = 0)), where masyn 
is the equilibrium value of the magnetization in a field h at the temperature T. Hence, 
measuring masyn "we can obtain the droplet theory estimate for the order parameter as: 

^droplet ^ ^ _ ^^^^ 
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Figure 5. Magnetization as a function of time for (3 = 1.25, L = 60 and h — 0.03. 
We have plotted in the inset the main plot with logarithmic scale in the abscissas, in 
order to make sure we had found a plateau. We have marked in the main figure the 
plateau of the magnetization with a horizontal line. 



We will shown in this section plots corresponding to (3 = 1.25 and L = 60. In order 
to obtain numerically masyn we have performed a very large in-field numerical simulation 
recording the value of the magnetization at the time t: m{t). The asymptotic value is 
simply masyn = m(cxo) (this observable shows really small dependence on L for the lattice 
sizes simulated in this paper). To avoid extrapolations we have continued the run until 
the magnetization shows a plateau (this means that the magnetization has reached its 
equilibrium value), and so we extract the value of masyn by computing the position of 
this plateau. For instance, we show in figure © the magnetization as a function of time 
for p = 1.25 and L = 60. 

By computing the asymptotic value of the magnetization for different temperatures, 
we obtain a reliable estimate for the order parameter in the droplet theory. In Table 
HI we report these values for the droplet theory estimates and, in addition, we write 
the values for the order parameter obtained in the first part of this paper, that we will 
denote in the rest of the paper as QeaW)- 

We recall that the values of q^AiP) reported in Table [H have small finite size effects 
(taking into account their error bars) as checked in figure Q- Moreover, we have found 
strong discrepancies between q'^^{(3) and ^ea^''^* small temperatures. 

We will describe in the rest of the paper our results for the violation of FDT out 
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p 




droplet 

Iea 


1.25 


0.6583(34) 


0.5573(13) 


1.00 


0.5071(31) 


0.3957(17) 


0.95 


0.3554(7) 


0.3404(21) 



Table 1. QEAiP) for L = 60 from C{t, t^) (obtained in the first part of the paper) and 
assuming droplet theory from mT/h. All the data showed in this table were obtained 
in a L = 60 lattice except for dynamical ^ea at /3 = 0.95 that was obtained simulating 
a L = 30 lattice. 

of equilibrium. 

In figure (jH)) we report the FD data out of equilibrium for one of the lowest 
temperature simulated. We have shown a vertical band which marks the our estimate 
of g^^, a straight line 1 — C to monitor the departure of this linear behavior and a 
horizontal band which marks masy^T/h (see figure (0))- In addition we have plotted 
data from three different waiting times. 

Figure © shows that our estimate for g^A^ matches very well in the plot and marks 
the region in which the FD data starts to depart from the linear behavior (for all the 
temperatures simulated). In figure (jZj) we have drawn a magnification of this region. 
In addition, in this figure one can see that the finite time effects in the building of the 
asymptotic curve are small. Practically the two biggest waiting times are compatible 
in the error (there is a factor ten in waiting time). With the state-of-the-art dedicated 
computed of the day it is impossible to simulate larger waiting times. We can conclude 
from this figure that we are unable to see dependence in waiting time for the two 
largest waiting times in the region in which they depart from the linear behavior. The 
dependence on the waiting time for larger times is smaller than our statistical errors. 
From our numerical data a droplet theory Fluctuation-dissipation asymptotic curve 
seems unlikely. 

6. CONCLUSIONS 

We have study numerically and out of equilibrium the three dimensional Ising spin glass 
with bimodal disorder. 

By computing the off equilibrium spin-spin correlation function we have been able 
to extract the order parameter of the phase transition. The study of the behavior of 
this order parameter with temperature permit us to compute the critical temperature 
and the associated critical exponent: both figures compare very well with previous 
numerical simulations and experiments. We have also discarded a XF-like scenario (we 
have found a non- vanishing order parameter in the low temperature region). We have 
also monitored the dependence of QEAiP) with the lattice size in the low temperature 
region for one (3. 

In the second part of the paper we have extracted the droplet prediction for the 




Figure 6. Fluctuation-dissipation curve out of equilibrium for one of the lowest 
temperature simulated (3 = 1.25, L = 60 and three waiting times. We have marked 
using three vertical lines the interval in which lies q'^^ for this [3 computed in the first 
part of the paper. In addition we have marked with three horizontal lines the value 
and the statistical error for rriasynT /h. Finally we have marked a vertical line with the 
droplet theory prediction for qea (left part of the plot) 



order parameter by computing the asymptotic value of the susceptibihty {mT/h). The 
droplet prediction compares (for all the /3's simulated) well with the order parameter 
computed in the first part of the paper for high temperature (of course, slightly below 
the critical temperature), but for lower temperatures the comparison is bad. 

Moreover the analysis (for larger waiting times) of the FD curves show a behavior 
that can be described in the RSB theory and points out that the droplet scenario seems 
unlikely (only a really small dependence on waiting time, outside of the precision of this 
work, could build a final FD curve compatible with the droplet theory). Moreover the 
point in which the numerical data depart from the linear behavior compares well with 
the estimate obtaining in the first part of this paper, supporting the RSB scenario. 
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Figure 7. Magnification of figure |S1 showing up the region in which the FD curves 
depart from the straight line 1 — C. We have marked using three vertical lines the 
interval in which lies computed in the first part of the paper. 
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